USE ElliDef
USE NodeInfoDef
USE metric
USE matrix_oper
USE Elliptic_Solvers
IMPLICIT NONE
TYPE (NodeInfo) :: Info
INTEGER N,error,idiag(10)
REAL (KIND=qPrec) :: dx,dy,dz
REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: p,derp,u,v,w,div,pl
REAL (KIND=qPrec), ALLOCATABLE, DIMENSION (:,:) :: ut,vt,wt
REAL (KIND=qPrec), ALLOCATABLE, DIMENSION (:,:,:,:) :: gg
! first thing allocate structures for sparse matrixes
INTEGER i,j,k,nx,ny,nz
REAL (KIND=qPrec) :: t0,t1







! populate Info structures
Info%mX=(/160,40,114/)

nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)

Info%nDim=3 ! 3D problem
! now define geometry of the problem
Info%dX=one
dx=one;dy=one;dz=one
Info%dmax=0.

! assume all Neumann conditions
Info%mthbc(1:6)=0


!now copy info to elli structure


ALLOCATE(Info%ElliInfo) ! allocate space for ellikit 
Info%ElliInfo%ndim=Info%ndim ! dimension of problem
Info%ElliInfo%mX=Info%mX ! size of grid
Info%ElliInfo%dX=Info%dX ! grid increments
Info%ElliInfo%mthbc=Info%mthbc ! definition of boundary conditions




error=CREATE_OPERATORS(Info%ElliInfo)


p=>Info%ElliInfo%p
CALL RANDOM_NUMBER(p)


div=>Info%ElliInfo%div

ALLOCATE(gg(0:nx+1,0:ny+1,0:nz+1,19))
CALL RANDOM_NUMBER(gg) ! to make sure the compiler is not doing dirty tricks

CALL CPU_TIME(t0)


CALL CPU_TIME(t0)
DO k=1,nz
   DO j=1,ny
      DO i=1,nx

         div(i,j,k)=gg(i,j,k,1)*p(i,j,k)+&
              gg(i,j,k,2)*p(i-1,j,k)+&
              gg(i,j,k,3)*p(i+1,j,k)+&
              gg(i,j,k,4)*p(i,j-1,k)+&
              gg(i,j,k,5)*p(i,j+1,k)+&
              gg(i,j,k,6)*p(i,j,k-1)+&
              gg(i,j,k,7)*p(i,j,k+1)+&
              gg(i,j,k,8)*p(i-1,j-1,k)+&
              gg(i,j,k,9)*p(i-1,j+1,k)+&
              gg(i,j,k,10)*p(i+1,j-1,k)+&
              gg(i,j,k,11)*p(i+1,j+1,k)+&
              gg(i,j,k,12)*p(i,j-1,k-1)+&
              gg(i,j,k,13)*p(i,j-1,k+1)+&
              gg(i,j,k,14)*p(i,j+1,k-1)+&
              gg(i,j,k,15)*p(i,j+1,k+1)+&
              gg(i,j,k,16)*p(i-1,j,k-1)+&
              gg(i,j,k,17)*p(i-1,j,k+1)+&
              gg(i,j,k,18)*p(i+1,j,k-1)+&
              gg(i,j,k,19)*p(i+1,j,k+1)
      END DO
   ENDDO
ENDDO
CALL CPU_TIME(t1)
   PRINT*,'time using do with g in array is',t1-t0
!
CALL CPU_TIME(t0)
idiag=(/0,1,nx-1,nx,nx+1,nx*ny-nx,nx*ny-1,nx*ny,nx*ny+1,nx*ny+nx/)
call mkl_ddiasymv('U', nx*ny*nz, gg(1,1,1,1), nx*ny*nz, idiag, 10, p(1,1,1), div(1,1,1))
CALL CPU_TIME(t1)
   PRINT*,'time using mkl_ddiasymv is',t1-t0

!
CALL CPU_TIME(t0)
div=>ADV_SYM(Info%ElliInfo%G(CSR_LAPL)%crsptr,p(1:nx,1:ny,1:nz),(/0,0,0/))
CALL CPU_TIME(t1)
   PRINT*,'time using ADV_SYM is',t1-t0

NULLIFY(p,div)
error=FREEUP_OPERATORS(Info%ElliInfo)
DEALLOCATE(Info%ElliInfo,gg)
STOP

END

